Using the sf and tmap packages to overlay a spatial object ontop of a leaflet map with an OS Maps API base map.


tmap

The tmap R package is providing a wrapper around the JavaScript web mapping library Leaflet.js. It allows the visualisation of geospatial data in R on an interactive leaflet map.

library(sf)
Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1; sf_use_s2() is TRUE
library(tmap)
Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio
The legacy packages maptools, rgdal, and rgeos, underpinning this package
will retire shortly. Please refer to R-spatial evolution reports on
https://r-spatial.org/r/2023/05/15/evolution4.html for details.
This package is now running under evolution status 0 

Create sf object from GeoPackage (GPKG)

# Create an sf data.frame object from the Greenspace file
osogs <- st_read('../../data/ordnance-survey/os-open-greenspace-gb.gpkg',
                 layer = 'greenspace_site')
Reading layer `greenspace_site' from data source 
  `/cloud/project/data/ordnance-survey/os-open-greenspace-gb.gpkg' using driver `GPKG'
Simple feature collection with 150149 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 9819.84 ymin: 8274.57 xmax: 655229.5 ymax: 1214133
Projected CRS: OSGB 1936 / British National Grid

Spatially subset data frame

Using coordinate-based indexer to spatially subset by bounding box (BBOX).

# Subset the data by intersection with BBOX

# Greater London bounding box
bbox <- st_bbox(c(xmin = 503568.1996, 
                  xmax = 561957.4962, 
                  ymin = 155850.7975, 
                  ymax = 200933.9026), 
                crs = st_crs(27700))

# Convert into a 'geometry'
bbox <- st_as_sfc(bbox)

# Spatial subsetting/filtering
osogs_filtered <- osogs[bbox, ]

OS Maps API ZXY resource

# Reproject to Web Mercator (EPSG: 3857)
osogs_filtered <- st_transform(osogs_filtered, crs = 3857)
# Basic plot using tmap
p <- tm_shape(osogs_filtered) +
  tm_fill(col = '#00cd6c') +
  tm_layout(frame = FALSE)
# OS Maps API layer
# Example uses Light Style in Web Mercator (EPSG:3857) projection
layer <- 'Light_3857'

# OS Data Hub project API key
key <- '7UTXMMWsGjjIzcBmLdAGnMO6WEAQi9Ng'

# Define the tile server parameters for the basemap
url <- paste0('https://api.os.uk/maps/raster/v1/zxy/', layer,
              '/{z}/{x}/{y}.png?key=', key)

Create custom interactive map

# Change plotting mode from static to interactive
tmap_mode('view')
tmap mode set to interactive viewing
# Combine the features and base map
final_plot <- p +
                tm_basemap(server = url)

final_plot
# Change back to the static plotting mode
tmap_mode('plot')
tmap mode set to plotting

Alternative method using leaflet

Demonstration using a different R package.

library(leaflet)
library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

Create custom map

# m is the map object all other layers will be added to.
m <- leaflet() %>%
  setView(lng = -0.105, lat = 51.507, zoom = 12)
m <- m %>% 
  addTiles(urlTemplate = url,
           attribution = paste0('Contains OS data © Crown copyright and database right ', 
                                format(Sys.Date(), "%Y")),
           group = 'OS Maps') %>%
  addPolygons(data = st_transform(osogs_filtered, 4326),
              color = '#00cd6c')

# Return map
m
LS0tCnRpdGxlOiAiMDYgQWRkIGFuIE92ZXJsYXkgdG8gYW4gSW50ZXJhY3RpdmUgTWFwIgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgpVc2luZyB0aGUgYHNmYCBhbmQgYHRtYXBgIHBhY2thZ2VzIHRvICBvdmVybGF5IGEgc3BhdGlhbCBvYmplY3Qgb250b3Agb2YgYQpsZWFmbGV0IG1hcCB3aXRoIGFuIE9TIE1hcHMgQVBJIGJhc2UgbWFwLgoKLS0tCgojIyBgdG1hcGAKClRoZSBbYHRtYXBgXShodHRwczovL3ItdG1hcC5naXRodWIuaW8vdG1hcC8pIFIgcGFja2FnZSBpcyBwcm92aWRpbmcgYSB3cmFwcGVyCmFyb3VuZCB0aGUgSmF2YVNjcmlwdCB3ZWIgbWFwcGluZyBsaWJyYXJ5IFtMZWFmbGV0LmpzXShodHRwczovL2xlYWZsZXRqcy5jb20vKS4KSXQgYWxsb3dzIHRoZSB2aXN1YWxpc2F0aW9uIG9mIGdlb3NwYXRpYWwgZGF0YSBpbiBSIG9uIGFuIGludGVyYWN0aXZlIGxlYWZsZXQKbWFwLgoKYGBge3J9CmxpYnJhcnkoc2YpCmxpYnJhcnkodG1hcCkKYGBgCgojIyBDcmVhdGUgYHNmYCBvYmplY3QgZnJvbSBHZW9QYWNrYWdlIChHUEtHKQoKYGBge3J9CiMgQ3JlYXRlIGFuIHNmIGRhdGEuZnJhbWUgb2JqZWN0IGZyb20gdGhlIEdyZWVuc3BhY2UgZmlsZQpvc29ncyA8LSBzdF9yZWFkKCcuLi8uLi9kYXRhL29yZG5hbmNlLXN1cnZleS9vcy1vcGVuLWdyZWVuc3BhY2UtZ2IuZ3BrZycsCiAgICAgICAgICAgICAgICAgbGF5ZXIgPSAnZ3JlZW5zcGFjZV9zaXRlJykKYGBgCgojIyBTcGF0aWFsbHkgc3Vic2V0IGRhdGEgZnJhbWUKClVzaW5nIGNvb3JkaW5hdGUtYmFzZWQgaW5kZXhlciB0byBzcGF0aWFsbHkgc3Vic2V0IGJ5IGJvdW5kaW5nIGJveCAoQkJPWCkuCgpgYGB7cn0KIyBTdWJzZXQgdGhlIGRhdGEgYnkgaW50ZXJzZWN0aW9uIHdpdGggQkJPWAoKIyBHcmVhdGVyIExvbmRvbiBib3VuZGluZyBib3gKYmJveCA8LSBzdF9iYm94KGMoeG1pbiA9IDUwMzU2OC4xOTk2LCAKICAgICAgICAgICAgICAgICAgeG1heCA9IDU2MTk1Ny40OTYyLCAKICAgICAgICAgICAgICAgICAgeW1pbiA9IDE1NTg1MC43OTc1LCAKICAgICAgICAgICAgICAgICAgeW1heCA9IDIwMDkzMy45MDI2KSwgCiAgICAgICAgICAgICAgICBjcnMgPSBzdF9jcnMoMjc3MDApKQoKIyBDb252ZXJ0IGludG8gYSAnZ2VvbWV0cnknCmJib3ggPC0gc3RfYXNfc2ZjKGJib3gpCgojIFNwYXRpYWwgc3Vic2V0dGluZy9maWx0ZXJpbmcKb3NvZ3NfZmlsdGVyZWQgPC0gb3NvZ3NbYmJveCwgXQpgYGAKCiMjIE9TIE1hcHMgQVBJIFpYWSByZXNvdXJjZQoKYGBge3J9CiMgUmVwcm9qZWN0IHRvIFdlYiBNZXJjYXRvciAoRVBTRzogMzg1NykKb3NvZ3NfZmlsdGVyZWQgPC0gc3RfdHJhbnNmb3JtKG9zb2dzX2ZpbHRlcmVkLCBjcnMgPSAzODU3KQpgYGAKCmBgYHtyfQojIEJhc2ljIHBsb3QgdXNpbmcgdG1hcApwIDwtIHRtX3NoYXBlKG9zb2dzX2ZpbHRlcmVkKSArCiAgdG1fZmlsbChjb2wgPSAnIzAwY2Q2YycpICsKICB0bV9sYXlvdXQoZnJhbWUgPSBGQUxTRSkKYGBgCgpgYGB7cn0KIyBPUyBNYXBzIEFQSSBsYXllcgojIEV4YW1wbGUgdXNlcyBMaWdodCBTdHlsZSBpbiBXZWIgTWVyY2F0b3IgKEVQU0c6Mzg1NykgcHJvamVjdGlvbgpsYXllciA8LSAnTGlnaHRfMzg1NycKCiMgT1MgRGF0YSBIdWIgcHJvamVjdCBBUEkga2V5CmtleSA8LSAnN1VUWE1NV3NHampJemNCbUxkQUduTU82V0VBUWk5TmcnCgojIERlZmluZSB0aGUgdGlsZSBzZXJ2ZXIgcGFyYW1ldGVycyBmb3IgdGhlIGJhc2VtYXAKdXJsIDwtIHBhc3RlMCgnaHR0cHM6Ly9hcGkub3MudWsvbWFwcy9yYXN0ZXIvdjEvenh5LycsIGxheWVyLAogICAgICAgICAgICAgICcve3p9L3t4fS97eX0ucG5nP2tleT0nLCBrZXkpCmBgYAoKIyMgQ3JlYXRlIGN1c3RvbSBpbnRlcmFjdGl2ZSBtYXAKCmBgYHtyLCBmaWcud2lkdGggPSAxMCwgZmlnLmhlaWdodCA9IDEwfQojIENoYW5nZSBwbG90dGluZyBtb2RlIGZyb20gc3RhdGljIHRvIGludGVyYWN0aXZlCnRtYXBfbW9kZSgndmlldycpCgojIENvbWJpbmUgdGhlIGZlYXR1cmVzIGFuZCBiYXNlIG1hcApmaW5hbF9wbG90IDwtIHAgKwogICAgICAgICAgICAgICAgdG1fYmFzZW1hcChzZXJ2ZXIgPSB1cmwpCgpmaW5hbF9wbG90CmBgYAoKYGBge3J9CiMgQ2hhbmdlIGJhY2sgdG8gdGhlIHN0YXRpYyBwbG90dGluZyBtb2RlCnRtYXBfbW9kZSgncGxvdCcpCmBgYAoKIyMgQWx0ZXJuYXRpdmUgbWV0aG9kIHVzaW5nIGBsZWFmbGV0YAoKRGVtb25zdHJhdGlvbiB1c2luZyBhIGRpZmZlcmVudCBgUmAgcGFja2FnZS4KCmBgYHtyfQpsaWJyYXJ5KGxlYWZsZXQpCmxpYnJhcnkoZHBseXIpCmBgYAoKIyMgQ3JlYXRlIGN1c3RvbSBtYXAKCmBgYHtyfQojIG0gaXMgdGhlIG1hcCBvYmplY3QgYWxsIG90aGVyIGxheWVycyB3aWxsIGJlIGFkZGVkIHRvLgptIDwtIGxlYWZsZXQoKSAlPiUKICBzZXRWaWV3KGxuZyA9IC0wLjEwNSwgbGF0ID0gNTEuNTA3LCB6b29tID0gMTIpCmBgYAoKYGBge3IsIGZpZy53aWR0aCA9IDEwLCBmaWcuaGVpZ2h0ID0gMTB9Cm0gPC0gbSAlPiUgCiAgYWRkVGlsZXModXJsVGVtcGxhdGUgPSB1cmwsCiAgICAgICAgICAgYXR0cmlidXRpb24gPSBwYXN0ZTAoJ0NvbnRhaW5zIE9TIGRhdGEgwqkgQ3Jvd24gY29weXJpZ2h0IGFuZCBkYXRhYmFzZSByaWdodCAnLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBmb3JtYXQoU3lzLkRhdGUoKSwgIiVZIikpLAogICAgICAgICAgIGdyb3VwID0gJ09TIE1hcHMnKSAlPiUKICBhZGRQb2x5Z29ucyhkYXRhID0gc3RfdHJhbnNmb3JtKG9zb2dzX2ZpbHRlcmVkLCA0MzI2KSwKICAgICAgICAgICAgICBjb2xvciA9ICcjMDBjZDZjJykKCiMgUmV0dXJuIG1hcAptCmBgYAo=